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1.  imODUCTION 


Errors  in  Inertial  navigation  can  be  subdivided  into  two 
groups,  the  system  related  errors,  and  the  modelling  errors*  The 
first  group  comprises  e.g.  measuring  noise,  calibration  and  alignment 
errors  while  the  second  group  contains  approximation,  linearization 
and  lumping  errors.  Traditionally,  advances  in  system  technology 
which  resulted  in  a  reduction  of  the  system  related  errors  have 
always  called  for  a  scrutiny  of  the  underlying  models.  For  a  long 
time,  one  of  the  approximations,  the  normal  gravity  model,  proved  to 
be  sufficient  for  airborne  Inertial  navigation.  However,  with  sensor 
accuracies  predicted  for  third  generation  inertial  systems  (Draper, 
1977),  this  will  not  be  true  anymore.  This  report  therefore  studies 
the  size  and  characteristical  behaviour  of  gravity  induced  position 
errors  in  airborne  inertial  navigation. 

In  order  to  relate  the  accuracy  of  a  given  gravity  field 
model  to  the  system  errors,  the  effect  of  model  approximations  on 
system  performance  must  be  studied.  This  topic  is  by  no  means  a  new 
one.  A  number  of  papers  in  the  late  sixties  and  early  seventies 
dealt  with  it  and  discussed  various  aspects  of  the  problem.  Since 
they  were  written  by  systems  engineers  and  were  published  in  journals 
seldom  read  by  geodesists,  they  had  almost  no  impact  on  geodesy. 

This  is  regrettable  for  two  reasons.  On  the  one  hand,  the  represen-- 
tation  of  a  velocity  transformed  gravity  field  adds  an  interesting  new 
dimension  to  the  usual  geodetic  approach.  On  the  other  hand,  geodetic 
results  on  the  covariance  structure  of  the  gravity  field  could  have 
been  readily  used  to  generate  models  consistent  with  our  present 
knowledge.  The  objective  of  this  report  is  therefore  twofold.  On 
the  one  hand,  it  will  present  the  basic  solution  approach  without 
assuming  a  background  in  systems  theory.  It  is  hoped  that  in  this 
way  some  interest  in  this  fascinating  area  will  be  stirred  among 
geodesists.  On  the  other  hand,  an  attempt  will  be  made  to  incorporate 
recent  geodetic  results  into  the  solution  and  to  point  out  areas  where 
further  research  is  needed. 


2.  TH9  DYNAMIC  ERROjEt  MDDRL 


Inertial  navlgatipn  is  based  on  Newton's  second  lav.  Thus, 
the  basic  nathematical  model  is  a  system  of  second-order  differential 
equations*  Such  a  system  can  always  be  transformed  into  a  system  of 
first-order  differential  equations*  The  derivation  of  the  first- 
order  model  from  the  basic  set  of  equations  is  discussed  in  detail  in 
Britting  (1971)  for  the  different  mechanizations  used  in  inertial 
navigation*  The  resulting  dynamic  system  usually  is  of  the  state 
space  variety,  i*e.  it  is  expressed  in  terms  of  a  time  dependent 
vector  x(t)  whose  cooqionents  at  any  specified  time  represent  the  state 
of  the  system*  The  dynamic  error  model  is  a  combination  of  a  state 
space  formulation  and  white  noise  processes.  Such  a  model  offers  a 
considerable  flexibility  in  applications  and  takes  into  account  that 
the  error  process  contains  both  deterministic  and  stochastic  elements* 
General  characteristics  of  the  dynamic  error  model  will  be  discussed 
in  this  sectipn  while  the  specific  formulation  of  the  Inertial  navi¬ 
gation  problem  will  be  the  topic  of  section  5* 

In  its  most  general  form  the  linear  dynamic  error  model  is 

given  by 

i(t)  -  F(t)x(t)  +  G(t)w(t)  +  L(t)u(t)  (2.1a) 

with  initial  conditions 

x(o)  -  c  (2*lb) 

and  control  measurements 

z(t)  •  H(t)x(t)  +  r(t)  (2.1c) 

where 

x(t)  ***  is  the  state  vector  of  the  system 

x(t)  *•*  is  the  time  derivative  of  x 

F(t)  • • •  is  the  dynamics  matrix 

G(t)w(t)  ..*  is  the  system  noise  (stochastic) 

L(t)u(t)  •*•  is  the  control  function  (deterministic) 
z(t)  • • •  is  the  measurement 


H(t)  •••  Is  the  design  matrix,  relating  measurement  and  state 
vector 

r(t)  •••  is  the  measurement  noise 
and  i^ere  all  underlined  quantities  represent  vectors  (lover  case)  or 
matrices  (upper  case).  In  these  equations  it  is  assumed  that  the 
matrices  F,  G,  H  are  completely  known,  i.e.  the  uncertainty  enters  via 
w(t),  r(t),  and  x(o).  The  variables  w(t)  and  r(t)  are  stochastic 
white  noise  processes  which  do  not  have  a  time  structure.  The  know¬ 
ledge  of  the  process  value  at  one  instant  of  time  does  not  provide 
knowledge  about  its  values  at  any  other  instant.  These  processes  are 
completely  defined  by  their  first  and  second  moments 

E{w(t)}  -  0  E{w(ti)w''’(t2)}  -  5(ti-t2)Q(t)  (2.2a) 

E{r(t)}  -  0  E{r(ti)r'^(t2)j  -  «(ti-t2)R(t)  (2.2b) 

where  E{  .  }  denotes  the  statistical  expectation  and  6(ti-t2)  is  the 
Dirac  delta  function.  This  definition  does  not  assume  that  the  process 
is  Gaussian  because  the  probability  densities  of  w(t)  and  r(t)  may  be 
specified  in  any  way.  The  initial  conditions  are  modelled  as  a  random 
vector 

E{x(o)}  =  0  E{x(o)x^(o)}  *  P(o)  (2.3) 

The  model  expressed  by  equations  (2.1)  to  (2.3)  has  the  important 
property  that  x(o) ,  w(t),  and  r(t)  are  not  related,  i.e. 

E{x(o)w^(t)}  -  0  (2.4a) 

E{x(o)r^(t)}  -  0  (2.4b) 

E{w(t)r^(t)}  *  0  (2.4c) 

The  deterministic  input  L(t)u(t)  is  usually  considered  as  known  which 
for  a  linear  model  like  (2.1)  means  that  its  effect  can  be  added  or 
subtracted.  It  will  therefore  be  disregarded  in  the  sequel.  A  further 
simplification  can  be  obtained  by  requiring  that  F,  G,  H,  Q,  and  R  are 
matrices  with  constant  coefficients.  In  this  case  a  time  invariant 
model  is  obtained  which  is  of  the  form 


x(t)  -  Jx(t)  +  Gw(t) 


(2.5a) 

(2.5b) 

(2.5c) 


with  x(o)  •  c 
and  z(t)  ■  ^(^)  +  £(t) 

It  should  be  noted  that  such  a  model  does  not  imply  that  x(t)  is  a 
stationary  stochastic  process.  The  advantage  of  model  (2.5)  versus 
model  (2.1)  is  that  statements  on  the  stability  of  the  system  are 
usually  simpler  and  that  it  is  sometimes  possible  to  find  an  analyti¬ 
cal  solution  to  the  homogeneous  system  of  differential  equations 

x(t)  -  Fx(t)  .  (2.6) 

Such  solutions  simplify  the  error  analysis  and  are  usually  the  only 

way  to  arrive  at  reliable  error  estimates  without  excessive  use  of 

computer  time.  Equation  (2.5)  will  therefore  be  considered  as  the 

basic  model  for  all  the  following  discussions. 

The  formal  solution  of  eqiiation  (2.5)  can  be  written  as 
t 

♦  (t,s)G^(s)ds  (2.7) 

0 

where  £(t,tQ)  is  called  the  fundamental  matrix  or  the  transition 
matrix.  It  satisfies  the  matrix  differential  equation 

i(t,to)  -  F(t)  £(t,to)  (2.8a) 

with  initial  conditions 

-  I  .  (2.8b) 

The  first  term  on  the  right-hand  side  of  equation  (2.7),  the  solution 
of  equation  (2.6),  is  called  the  complementary  function  or  the  tran¬ 
sient  solution.  It  gives  the  solution  of  the  system  of  differential 
equations  in  absence  of  an  input  or  forcing  function  Gw(t).  The  second 
term  is  called  the  particular  integral  or  the  steady  state  solution. 

It  adds  to  the  solution  the  term  originating  from  the  functions  Gw(t) 
in  equation  (2. 5a)* 

For  a  constant  coefficient  matrix  F  the  transition  matrix  is 
given  by  the  matrix  exponential 


x(t)  -  «(t,tg)c  + 


/ 


-  5  - 


«(t,to)  -  (2.9a) 

or  using  the  well-known  series  expansion  of  the  exponential  by 
.  (t-to)^ 

♦  (t.to)  -  I  - 7- -  (2.9b) 

1-0  i  ! 


Using  the  Cayley-Hamllton  theorem  the  Infinite  series  expansion  can 
be  replaced  by  a  finite  series  of  order  N  where  N  Is  the  rank  of  F. 

A  discussion  of  this  approach  and  its  usefulness  in  applications  is 
given  in  Blermann  (1971).  However,  it  is  not  always  the  most  direct 
way  to  a  solution.  Often  the  use  of  the  Laplace  transform  is  simpler* 
The  Laplace  transform  L  of  a  function  x(t)  is  defined  as 


f(s) 


-  L{x(t)} 

00 

■  e~^^x(t)dt 
0 


(2.10a) 

(2.10b) 


where  t>0  and  the  Integral  converges  for  some  value  of  s.  Similarly, 
for  a  vector  x(t)  of  functions  x^(t) 

f(s)  -  L{x(t)} 

00 

=  f  e“®‘^x(t)dt  (2.11b) 

0 


where  the  operation  (2.10b)  is  now  performed  on  each  element  of  x(t). 
Taking  the  Laplace  transform  of  equation  (2.6)  results  in 


sf(s)  -  x(0)  -  F  f(s)  (2.12a) 

and  using  (2.5b)  we  obtain 

(si  -  1)  f(s)  -  c  .  (2.12b) 

Thus,  the  system  of  differential  equations  (2.6)  has  been  converted 
into  a  system  of  algebraic  equations.  If  the  inverse  of  (si  -  F) 
exists  we  can  write 

f(s)  -  (si  -  P®c 


where  the  superscript  g  denotes  some  suitably  defined  inverse.  In 
the  following »  we  will  only  use  the  Cayley  Inverse,  thus  obtaining 

f(s)  -  (si  -  F)”^c  .  (2.12c) 

Hochstadt  (1975)  discusses  In  detail  the  existence  of  the  Cayley  in*^ 
verse  for  the  above  case. 

The  Inverse  Laplace  transform  L  ^  reverses  the  operation  (2.11), 

l.e. 

x(t)  =»  L  ^{f(s)}  (2.13) 


where  the  same  symbols  as  in  equation  (2.11)  have  been  used.  The 
functions  x(t)  and  f(s)  are  often  called  a  Laplace  transform  pair. 

Thus,  the  solutions  to  the  algebraic  equations  (2.12c)  can  be  inverse 
transformed  to  provide  solutions  to  the  original  differential  equa¬ 
tions  (2.6).  This  elegant  and  powerful  technique  will  be  frequently 
used  In  the  following  and  therefore  a  simple  example  will  be  given  here 
to  demonstrate  the  salient  features  of  the  technique. 

Consider  the  system  of  differential  equations 

6v  =  -K6v  -  u)^6p 
(Sp  =  5v 


which  Is 

of 

type  (2.6) 

with 

_ 

r*  -  “1 

6v 

N 

3 

1 

1 

X  * 

.  F  - 

and  c  = 

i 

6p 

1  0 

We  thus  have 


S  +  K  0)2 

(si  -  F)  » 

-1  s 


and  Its  Inverse 
(si  -  - 


1 

s  (s  +  K)  +  u)^ 


s 

1 


s  +  K 


(2.14) 


The  Individual  components  of  f(s)  are  therefore 


we  obtain 


L  -  x’j(t)  +  f2i(0)6(t) 


or 


I  s(s  +  K)  +  0)^  j  D 


e  D  cos  D  t  -  y  sin  D  t}  + 


(2.15a) 


x^j(t)  . 


Finally  £22(8)  bd  split  Into  the  two  components 
^  ^  8  K 

f  (s)  - -  +  - - 

8(8 +  K)  +  0)2  8(8 +  K)  +  0)2 

which  can  be  Laplace  transformed  using  equations  (2.15a)  and  (2.15c) 


Thus  the  solution  of  equation  (2.14)  can  be  written  as 
x(t) 


6v 

3 

CjXjj(t)  -  CjXjjCt) 

_<5P_ 

CiX^iCt)  + 

(2.16) 


where  the  terms  3t^^(t)  to  X22(t)  are  given  by  formulas  (2.15a)  to 
(2.15d). 

This  simple  example  shows  already  that  It  Is  important  to  find 
the  inverse  Laplace  transform  of  the  determinant  of  the  Inverse  matrix. 
The  Inverse  Laplace  transform  of  most  elements  of  (si  -  F) can  then 
be  obtained  by  rather  elementary  manipulations.  The  example  does  not 
show  that  the  analytical  Inversion  of  (si  -  F)  Is  a  major  problem.  If 
the  size  of  the  matrix  exceeds  10,  It  is  often  Impossible  to  come  up 
with  manageable  exptesslons. 


-  9  - 


The  basic  model  (2, 5}  Is  not  always  sufficient  to  cover  all 
cases  encountered  in  practical  problems.  Two  important  extensions 
are  the  modelling  of  coloured  noise  and  nonlinearities  in  the  basic 
model.  Time  correlated  noise  processes  will  be  tised  to  model  the 
effect  of  the  anomalous  gravity  field.  The  usual  approach  is  by 
state  vector  augmentation.  The  coloured  noise  w^  is  modelled  as  the 
output  of  a  linear  system 

+  lia  (2.17) 


\diere  Is  again  white  noise.  This  equation  is  added  to  eqxiatlon 
(2.5)  resulting  In 


(2.18) 


where  w  In  equation  (2.5)  Is  defined  by 


W  **  W  w 


and  w^  is  white  noise.  Because  of  equation  (2.18)  this  is  called  state 
vector  augmentation.  The  assumption  is  that  the  correlated  noise  pro¬ 
cess  can  indeed  be  modelled  by  (2.17).  In  theory,  this  is  not  always 
possible.  In  practice,  it  can  usually  be  done  with  sufficient  accuracy. 
In  engineering  applications  the  following  F^-matrices  are  often  used: 

F  »  -S  , 
w  * 

modelling  a  first-order  Markov  process,  and 
Fw  -  0  . 

modelling  a  bias  term.  A  more  detailed  discussion  of  suitable  models 
for  the  anomalous  gravity  field  will  be  given  in  section  4. 

Nonlinearities  in  the  basic  model  will  play  no  role  in  the 
following  and  reference  is  therefore  made  to  standard  textbooks  as 
Gelb  (1974)  and  Eykhoff  (1974). 
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3.  ERBOR  PROPAGATION  AND  STEADY  STATE  ERRORS 

There  are  two  major  practical  advantages  of  using  the  dynamic 
model  (2.5) •  First*  the  real-time  implementation  of  these  formulas 
including  update  measurements  presents  no  problem.  Kalman  filtering 
and  its  alternatives  have  been  investigated  in  great  detail  and  their 
numerical  characteristics  are  well-known.  Second*  a  differential 
equation  for  the  error  propagation  can  be  derived  from  which  an  ex¬ 
pression  for  the  steady  state  error  of  individual  states  can  be 
obtained.  The  derivation  of  the  basic  error  eqtiation,  the  linear 
variance  equation*  and  its  solution  for  the  steady  state  errors  will 
be  given  in  this  section.  The  application  of  this  technique  to 
derive  analytical  expressions  for  the  position  errors  in  airborne  Iner' 
tial  navigation  will  be  given  in  section  5. 

The  statistical  description  of  the  total  error  of  a  dynamical 
system  at  a  certain  time  t  is  given  by  the  covariance  matrix  P(t) 

P(t)  -  E{x(t)  /(t)}  .  (3.1) 

The  change  of  P(t)  with  time  models  the  error  propagation.  We  will 
thus  consider  the  function  P(t»  t-e) 

P(t,  t-e)  »  E{x(t)  x^(t-e)}  (3.2) 

where  e  is  a  small  positive  number*  and  then  use  limit  arguments  to 
arrive  at  the  result.  Differentiation  of  equation  (3.2)  with  respect 
to  time  yields 

P(t*  t-e)  »  E{x(t)x^(t-e)  +  x(t)x^(t-€)}  . 

Using  equation  (2.5a)  we  obtain 

P(t*  t-c)  *  E{  [Fx(t)  +  ^(t)]x^(t-e)  + 

x(t)I^'^(t-e)F^  +  w^(t-e)G^]} 
or 

P(t,  t-e)  -  E{FX(t)x^(t-e)  +  Gw(t)x^(t-e)  + 

x(t)x^(t-c)F^  +  x(t)w^(t-e)G^}  . 


(3.3) 
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Because  of  equation  (2.4a)  the  statistical  expectation  of  the  second 
term  in  equation  (3.3)  must  be  zero.  The  first  and  third  term  are 
easy  to  evaluate.  The  last  term  is  rewritten  using  equation  (2.7) 
for  x(t)  and  equation  (2.5b)  for  c 

E{x(t)w^(t-e)G^}  •  E{*(t,o)x(0)w^(t-c)  + 


n 


*(t,s)Gw(s)w^(t-e)G^ds} 


Again,  the  evaluation  of  the  first  term  results  in  zero  because  of 
equation  (2. 4a)  while  the  second  term  can  be  written  as 

E{x(t)w^(t-€)G^}  -  £(t,  t-e)  GQ  (t-e)G^ 

using  equation  (2. 2a)  •  Thus  we  obtain 

P(t,  t-e)  -  E{^(t)x'^(t-e)  +  x(t)x’^(t-e)F'’^  + 


£(t,  t-e)  Gg  (t-e)G’^} 
Using  e  ->■  0,  equation  (3.1)  and 
«(t,t)  -  I 


(3.4) 


we  get 

P(t)  -  FP(t)+P(t)F^  +  Gg(t)G^  (3.5) 

which  is  the  required  expression.  Its  name,  linear  variance  equation^ 
is  due  to  the  fact  that  it  is  linear  in  P.  If  update  measurements 
(2.5c)  are  made,  a  similar  derivation  leads  to  the  matrix  equation 

P(t)  -  FP(t)  +  P(t)  F^  +  Gg(t)G^  -  P(t)  HR“^HV(t)  (3.6) 

which  obviously  is  nonlinear  in  P.  This  equation  is  known  as  the 
matrix  Riccati  equation  and  has  been  thoroughly  investigated  in  the 
mathematical  literature  in  connection  with  variational  problems.  Its 
solution  is  generally  obtained  by  transforming  the  nonlinear  equation 
into  a  set  of  two  linear  equations  of  the  form 


X(t) 

■-F^ 

T  -1 
HR  H 

X(t) 

Y(t) 

GQG^ 

F 

T(t) 
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with  initial  values 


x(to)  -  I 

The  matrices  Y  and  X  are  related  via 


(3.7b) 


Y(t)  -  P(t)  X(t)  (3.8) 

from  \dilch  P(t)  can  be  obtained 

P(t)  -  Y(t)  X'Vt)  .  (3.9) 

These  transformations  can  be  verified  by  differentiating  equation 
(3.9)  and  transforming  It  Into  the  second  equation  (3.7a)  by  use  of 
equation  (3.6)  and  the  first  row  In  (3.7a).  For  details,  see 
Llebelt  (1967).  X(t)  Is  the  transition  matrix  of  the  differential 
equation 

X(t)  -  {-F^  +  (t)}  X(t)  (3.10) 


which  Is  the  adjoint  of  the  differential  equation  (3.6).  Since  the 
Inverse  of  a  transition  matrix  always  exists,  the  validity  of  equa¬ 
tion  (3.9)  Is  secured. 

To  obtain  the  matrices  X  and  Y,  we  write  the  formal  solution 
of  the  system  (3.7a)  by  way  of  a  transition  matrix  ♦(t,tQ)  as 


X(t)“ 

Y.(t) 


and  using  equation  (3.7b)  and  partitioning 


X(t) 

Y(t) 


£ii(t,to)  £i2(t,to) 
i2l(t,to)  *22(t,to) 


where  again  p 


4(t.tg) 


-P‘ 


992 


T  -1 
HR  H 


I 

-1 

P(to) 

0^ 

(3.11) 


(3.12) 

1,2)  can  Immediately  be 


The  dimensions  of  the  submatrices  (l,j 
obtained  from  equation  (3.12).  Using  equation  (3.9)  and  (3.11)  we 
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finally  get  the  solution 

P(t)  -  {i2l(t,to)+*22(t.t5)P(to)}{*ii(t,tQ)+«i2(t,tQ)P(tj,)}  (3.13) 


which  gives  the  covariance  matrix  P  at  an  arbitrary  time  t  in  terms 
of  its  initial  values  P(tQ)  and  the  transition  matrices  of  the  system 
(3e7a)e  This  solution  goes  back  to  Reid  (1946).  A  proof  and  general 
discussion  is  given  in  Levin  (1959).  Its  first  application  to  optimal 
estimation  problems  can  be  found  in  Kalman  and  Bucy  (1961). 

The  practical  use  of  equation  (3«13)  is  somewhat  limited 
because  the  transition  matrices  are  seldom  given  in  analytical  form. 
Usually,  they  are  computed  by  numerical  integration  and  in  this  case 
recurrence  relations  are  much  more  efficient  to  use.  In  our  case, 
however,  equation  (3.13)  can  be  sufficiently  simplified  to  obtain  use- 
ftxl  analytical  expressions.  The  first  simplification  is  due  to  the 
fact  that  we  do  not  have  to  solve  the  general  equation  (3.6)  but  the 
special  case  (3.5).  Noting  that 
T  -1 

HR  H  -  0 


and  using  the  Laplace  transform  technique  on  equation  (3.12),  we 
obtain 


♦(t.tg)  -  L"^  < 


and  after  performing  the  matrix  inversion 


T 

-1 

< 

sI+F  0 

> 

-G  g  si  -  F 

1 

4 

*<t,tn)  -  L 


-1 


(sI  +  F^)"^ 


(si  -  F)"^  G  (sl+ j  (sI-F) 


-1 


(3.14) 


and  for  the  individual  submatrices 

*l^(t,tQ)  -  l"^{(sI+F^)"^} 


(3.15a) 

(3.15b) 


\ 
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JjjCt.tfl)  -  L"^{(8l-F)"^GgG'^(8l  +  F’^)"^}  (3.15c) 

£22(t,to)  -  L"^{(8l-n‘^>  .  (3.15d) 

Furthermore,  we  are  only  interested  in  the  steady  state  errors,  i.e. 
in  that  part  of  P(t)  which  is  independent  of  P(tQ)«  Thus,  observing 
equation  (3.15b)  and  the  steady  state  requirement,  equation  (3.13) 
simplifies  to 

Pgg(t)  »  •  (3.16) 

where  the  subscripts  'ss'  denote  the  steady  state  solution.  The  matrix 
identity 

(sI  +  f’^)"^  -  {sI  +  F)”^’^  (3.17) 

can  be  api^lied  to  show  that 

♦IjCt.tp)  -  *H(t.to)  (3.18) 

and  is  also  useful  in  simplifying  the  algebraic  operations  in  (3.15c). 
The  diagonal  terms  of  P  will  be  used  in  section  5  to  derive  analy- 

“SS 

tical  expressions  for  the  steady  state  position  error. 


1 


I 


I 
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4.  COVARIANCE  MODELS  FOR  THE  ANOMALOUS  GRAVITY  FIELD 

In  order  to  apply  the  error  analysis  developed  in  the  last 
section,  the  effects  of  the  anomalous  gravity  field  have  to  be 
modelled  into  a  form  compatible  with  equation  (2.5).  This  requires, 
first  of  all,  that  the  characteristical  features  of  the  field  can  be 
described  in  statistical  language,  e.g.  in  terms  of  random  functions 
and  their  associated  covariance  functions.  The  justification  of  such 
an  approach  is  given  in  Moritz  (1980)  and  within  the  limits  indicated 
there,  it  is  applicable  to  the  following  analysis.  The  second 
requirement  is  the  transformation  of  the  position  dependent  covari¬ 
ance  function  into  a  time  dependent  covariance  function  by  way  of  the 
velocity.  As  long  as  the  velocity  is  constant  and  the  basic  co- 
variance  function  isotropic,  this  transformation  corresponds  to  a 
simple  "time-scaling"  of  distances  between  points.  Since  a  constant 
aircraft  velocity  is  a  reasonable  assumption  in  our  case,  and  since 
all  gravity  anomaly  covariance  functions  in  practical  use  are  iso¬ 
tropic,  this  transformation  does  not  present  a  problem.  More  impor¬ 
tant  for  the  model  selection  are  two  other  requirements.  One  is  that 
the  model  is  self-consistent,  the  other  that  it  can  be  expressed  by  a 
set  of  differential  eqxiations  of  the  form  (2.17).  Self-consistency 
means  that  the  mathematical  relationships  which  are  known  to  exist 
between  the  anomalous  potential  and  linear  functionals  of  it  have  to 
be  used  when  deriving  the  covariance  functions  of  these  quantities. 

In  other  words,  the  covariance  functions  of  the  potential  and  the 
covariance  functions  of  its  derivatives  are  not  independent.  In  the 
geodetic  literature  self-consistency  with  respect  to  the  covariance 
functions  is  usually  discussed  in  the  context  of  covariance  propaga¬ 
tion.  The  other  requirement  is  obviously  important  whenever  an  error 
model  of  type  (2.5)  is  used.  It  has  found  considerable  attention  in 
the  navigation  literature  but  has  so  far  been  neglected  by  geodesists. 

For  our  analysis  both  aspects  are  equally  important. 

When  comparing  the  resulting  covariance  models,  two  major 
groups  can  be  distinguished.  They  will  be  labeled  'geodetic'  and 

•1 

1 

j 
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'navigational'  to  indicate  their  predominant  user  groups.  The 
geodetic  models  are  based  on  spherical  approximation  and  are  self- 
consistent  and  global  in  character.  Upward  continuation  of  these 
functions  is  no  problem  and  the  derivation  of  local  models  from  the 
global  function  is  relatively  simple.  The  important  parameters  of 
these  models  are  derived  from  different  data  groups,  specifically 
spherical  harmonic  coefficients  obtained  from  satellite  observations, 
gravity  anomalies  and  second-order  gradients.  Their  major  disadvan¬ 
tage  for  an  error  analysis  of  this  type  is  that  there  is  no  simple 
way  to  express  them  in  the  form  (2.17).  A  detailed  discussion  of  the 
mathematical  structure  of  these  functions  is  contained  in  Tscherning 
(1976)  and  Moritz  (1976).  The  fitting  of  a  global  function  is 
described  in  Tscherning  and  Rapp  (1974)  and  Jekeli  (1978)  while  local 
aspects  are  treated  in  Schwarz  and  Lachapelle  (1980)  and  Lachapelle 
and  Schwarz  (1980) .  Efficient  numerical  methods  to  compute  the  co- 
variance  function  are  treated  in  Sunkel  (1979,  1980).  Some  modelling 
aspects  of  nonhomogeneous  global  covariance  functions  are  discussed 
in  Rummel  and  Schwarz  (1977) .  The  navigational  models  are  based  on 
a  flat  earth  approximation  and  are  local  in  character.  Their  major 
advantage  is  that  they  can  be  cast  into  the  form  (2.17),  i.e*  they 
can  be  represented  as  the  output  of  a  linear  system  driven  by  white 
noise.  Not  all  of  them  are  self-consistent  and  most  of  them  cannot 
be  upward  continued  in  a  simple  way.  Usually  the  upward  continuation 
integral  has  to  be  evaluated.  The  extension  of  the  local  function  to 
a  global  function  is  in  general  not  possible.  The  important  parameters 
of  these  functions  are  usually  determined  from  one  data  group  only, 
specifically  gravity  anomalies  or  deflections  of  the  vertical. 

Levine  and  Gelb  (1969)  were  the  first  to  model  deflections  of  the 
vertical  in  this  way.  Subsequent  contributions  are  due  to  Shaw  et  al 
(1969),  Vyskosil  (1970),  Grafarend  (1971),  Kasper  (1971)  but  it  was 
especially  Jordan  (1972,  1973,  1981  et  al)  who  in  a  series  of  papers 
discussed  both  the  geodetic  and  navigational  aspects  of  different 
models  and  synthesized  the  work  of  his  predecessors. 
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This  brief  discussion  shows  already  that  for  the  airborne 
case  a  combination  of  the  two  approaches  Is  needed*  On  the  one  hand, 
the  geodetic  model  provides  a  simple  procedure  to  compare  different 
approximations  of  the  gravity  field  In  a  consistent  manner  and  to 
evaluate  the  effect  of  upward  continuation*  On  the  other  hand,  only 
the  navigational  model  can  be  used  In  a  state  space  formulation  which 
forms  the  basis  for  the  dynamic  error  model*  All  covariance  functions 
will  therefore  In  the  sequel  be  based  on  a  consistent  geodetic  model 
and  only  In  the  last  step  the  transition  to  the  navigational  model 
will  be  made*  This  last  step  caxinot  be' done  rigorously  but  will 
always  Involve  an  approximation.  To  keep  approximation  errors  small, 
the  fitting  of  the  navigational  model  will  be  done  by  way  of  the 
essential  parameters  defined  by  Moritz  (1976)*  Agreement  of  these 
parameters  for  two  models  secures  agreement  of  the  error  estimates 
derived  from  them*  The  remainder  of  this  section  will  review  the 
formulation  of  the  geodetic  model,  define  the  essential  parameters, 
and  discuss  some  simple  navigational  models  which  can  be  used  as 
approximations « 

The  rotation  Invariant,  spatial  covariance  function  of  the 
anomalous  potential  Is  given  by  the  general  expression 


K(P,Q) 

where 

P,Q 

rp.rq 

Pn 

'f'PQ 


AEkn 


n+1 

p„(cos 


(4.1) 


•  are  two  points  In  space  outside  the  sphere  with  radius  R, 
and  with  origin  in  the  centre  of  gravity  of  the  earth 

•  are  positive  coefficients,  l.e*  ^  0 

•  are  radial  distances  from  the  centre  of  the  sphere 

•  are  the  Legendre  polynomials 

•  Is  the  spherical  distance  between  P  and  Q. 


R(P,Q)  Is  a  function  of  positive  type,  symmetric  and  harmonic  In  the 


space  outside  the  sphere  R  and  regtilar  at  Infinity;  for  a  discussion, 
see  Krarup  (1969)*  The  coefficients  can  be  defined  rather 


-  18  - 


arbitrarily  as  long  as  their  positivity  Is  secured.  However,  the 
number  of  choices  Is  considerably  reduced  when  asking  for  coeffi¬ 
cients  which  result  In  a  simple  closed  expression  for  equation  (4.1). 
This  requirement  Is  necessary  for  any  practical  application  and  has 
led  to  coefficients  of  the  general  type 

k  -  n  (ti  +  b.)  ^ 

“  j-o  J 

where  H  denotes  the  multiplication  of  the  1  factors  on  the  right-hand 
side.  Well-known  examples  are 


or  1 


bo  =  -2 


resulting  In  coefficients 


(n-1)  (n-2) 


(4.3a) 


(n-l)(n-2)(n+24) 


(4.3b) 


By  expressing  the  as  a  sum  of  partial  fractions  the  Infinite  sum 
In  equation  (4.1)  can  be  expressed  as  a  sum  of  closed  analytical  ex¬ 
pressions  resulting  from  the  summation  of  the  corresponding  1 
Infinite  series.  As  an  example  the  covariance  function  of  the  model 
(4.3b)  can  be  written  as 


K(P,Q) 


tdiere 


125  -  26  [F.,  -  s3p^(t),  + 

sit  s3p2(t)  ^ 

■  25  26  ^ 


s(M  +  St  In  ^) 

i  s{M(l  +  3st)}  +  s^p, (t)  In  4  + 


'(1  -  t2) 


(4.4) 


(4.5a) 


(4.5b) 
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and  can  be  computed  from  the  recursion  formula 

Vl  "  irf;  +  (2k-l)tF^  -  (k-l)s-4j^_j}  .  (4.5c) 

The  auxiliary  quantities  L,M,N  are  defined  by 

L  -  (1  -  2st  + 

M  -  1  -  L  -  St 
N  -  1  +  L  -  St 
where 

t  -  cos 
r2 

s  ■  - 

r  r 

P  Q 

and  A  and  P2(t)  have  been  defined  earlier.  For  a  detailed  derivation 
reference  is  made  to  Tscheming  and  Rapp  (1974).  Expressions  of  the 
form  (4.5a)  to  (4.5c)  are  easy  to  evaluate  on  a  computer  and  are 
simple  enough  to  allow  the  formation  of  all  derivatives  necessary  for 
covariance  propagation.  Thus,  the  basic  model  given  by  equations 
(4.1)  and  (4.2)  is  used  to  derive  a  self-consistent  set  of  covariance 
functions  for  functionals  of  the  anomalous  potential.  These  deriva¬ 
tions  are  given  in  Tscheming  and  Rapp  (1974)  and  Tscheming  (1976) 
for  a  number  of  different  models. 

To  obtain  a  satisfactory  fit  of  empirically  determined  data  to 
equation  (4.1),  the  following  variables  can  be  changed: 

•  the  size  of  A 
.  the  radius  of  the  sphere  R 

.  the  coefficients  according  to  one  of  the  models  (4.2). 

The  effect  of  a  change  in  the  first  variable  is  simply  that  of  a 
scaling  factor  for  all  coefficients  and  thus  for  the  total  covariance 
function.  It  is  directly  related  to  the  variance  K(0,0)  or  the 
derived  variance  of  the  gravity  anomaly  C(0,0)  which  can  be  easily 
determined  from  given  point  gravity  anomalies.  A  change  of  R  affects 
the  spectrum  of  K(P,Q)  in  a  non-uniform  way.  For  instance,  a  decrease 
of  R  works  similar  to  a  low  pass  filter  by  strongly  damping  the 
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higher  frequencies.  Similarly,  when  going  from  model  (4.3a)  to 
model  (4.3b)  a  damping  of  the  high  degree  terms  will  be  the  result. 

It  Is  therefore  Important  to  ensure  that  a  covariance  function  of 
type  (4.1)  fits  well  to  data  representing  different  parts  of  the 
spectrum.  Typically,  geoidal  undulations  or  geopotential  models 
from  satellite  observations  contain  mainly  low  frequency  information 
while  torsion  balance  measurements  or  gradiometer  data  represent  the 
high  frequency  range.  Gravity  anomalies  and  deflections  of  the 
vertical  take  an  intermediate  position  containing  mostly  information 
on  the  medium  frequency  range.  At  present  it  appears  that  a  model 
of  type  (4.3b)  or  a  combination  of  this  model  with  one  of  a  simpler 
structure  gives  the  best  fit  to  all  available  data.  It  was  therefore 
decided  to  use  model  (4.3b)  for  the  following  computations.  Once  the 
basic  structure  of  the  model  is  given,  we  have  to  specify  what  will 
be  understood  by  goodness  of  fit.  For  instance,  is  it  important  to 
approximate  an  empirical  covariance  function  as  closely  as  possible 
everywhere  or  are  certain  regions  especially  important?  Or,  should 
empirical  degree  variances  computed  from  satellite  solutions  replace 
the  low  order  in  eqiiation  (4.1)?  Questions  of  this  type  have 

been  answered  by  Moritz  (1976,  1977)  in  a  discussion  of  the  essential 
parameters  of  planar  covariance  functions.  These  functions  can  be 
considered  as  planar  equivalents  of  the  covariance  functions  (4.1) 
and  results  obtained  for  them  carry  over  to  the  spherical  case  with 
only  small  modifications. 

Moritz  (1976)  defines  three  essential  parameters  for  co- 
variance  functions  of  this  type:  the  variance  Kq,  the  correlation 
length  t,  and  the  curvature  parameter  x  •  The  variance  Kq  is  the 
value  of  the  covariance  function  K(iji)  for  the  argument  *  0,  i.e. 

Kq  -  K(0)  (4.6a) 

The  correlation  length  is  the  value  of  the  argument  for  which 

K(?)  -  J  Kg  * 


I 
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The  curvature  parameter  x  is  related  to  the  curvature  k  of  K(i^)  at 
■  0  by 


X  -  K  •  C^Kn 


where  k  Is  given  In  the  usual  way  by 
K" 


(1  + 


and 


K* 


with 


3p 

R\|i  . 


f"  - 

3«2 


(4.6c) 

(4.7) 


Usually,  the  fit  Is  made  to  empirical  values  of  the  gravity  anomaly 
covariance  function  C(\|;)  because  large  data  samples  are  available. 

In  this  case,  the  curvature  k  is  equal  to  Gq  ,  the  variance  of  the 
horizontal  gradients  of  Ag  .  Thus  we  have 

X  -  /  Cq  ;  (4.8) 

for. a  derivation  see  Moritz  (1976).  In  this  case  all  three  essential 
parameters  can  be  determined  from  actxial  measurements.  is  the 
global  variance  of  the  gravity  anomalies;  ?  is  the  correlation  length 
of  the  empirical  covariance  function  C(^)  which  can  be  determined 
graphically;  Gq  can  be  computed  either  from  dense  point  gravity  anoma¬ 
lies  or  from  torsion  balance  or  gradiometer  measurements.  For  a 
more  detailed  discussion,  see  Schwarz  and  Lachapelle  (1980) * 

Two  ccrariance  functions  which  have  the  same  essential  para¬ 
meters  are  not  necessarily  equal  everywhere.  However,  their  general 
behaviour  with  respect  to  interpolation  will  be  similar.  This  is 
obviously  a  very  important  property  for  the  analysis  in  this  report. 

It  can  also  be  shown  that  the  shape  of  the  covariance  curve  at  dis¬ 
tances  larger  than  about  1.5 C  does  not  influence  the  interpolation 
very  much.  This  is  quite  important  for  any  empirical  covariance 
fitting. 

Besides  relating  the  covariance  model  (4.1)  to  the  actual 
data,  the  set  of  essential  parameters  can  also  be  used  to  get  a  good 


fit  between  the  geodetic  and  the  navigational  models.  This  will  be 
the  main  application  in  the  following.  Covariance  functions  for  the 
deflections  of  the  vertical  will  be  derived  from  the  best  available 
model  (4.1).  Upward  continuation  in  this  model  requires  only  a  change 
of  tp  and  r^  •  Improved  gravity  models  can  be  simulated  by  subtract¬ 
ing  a  certain  number  of  low  degree  coefficients  from  the  series.  By 
basing  all  computations  on  one  fundamental  model  of  type  (4.1),  self- 
consistency  between  different  functionals  of  T  as  well  as  consistency 
between  the  different  cases  is  secured.  By  defining  each  covariance 
function  in  terms  of  its  essential  parameters,  self-consistency  of 
the  approximating  navigational  model  is  not  required  anymore.  Emphasis 
will  therefore  be  on  selecting  navigational  models  which  give  a  best 
fit  in  terms  of  essential  parameters  for  the  individual  covariance 
functions. 

The  navigational  covariance  models  most  frequently  used  are 
low-order  Gauss-Harkov  models  whose  covariance  ftmctions  are  of  the 
general  form  (Gelb^  1974) 


K(t) 


e 


-Bnl- 


n-1 

Z 

k-0 


r(n)(2B^|T|) 


n-k-1 


+  m^ 


(4.9) 


(2n-2)  !  k  I  r(n-k) 


where 


m 


r(n) 


is  the  variance  of  the  process 

is  a  parameter  related  to  the  correlation  length 

is  the  mean 

is  the  Gamma  function 


n  ...  is  the  order  of  the  Gauss-Markov  process 

and  |t|  is  the  independent  variable,  l.e.  distance  in  geodetic  appli¬ 
cations  and  time  in  navigation  applications.  Implicit  in  this  defi¬ 
nition  is  that  the  process  is  stationary  and  its  probability  density 
function  is  Gaussian.  Processes  of  this  type  are  defined  by  a 
stochastic  differential  equation  of  order  n  or  equivalently  by  a 
set  of  n  first-order  stochastic  differential  equations.  They  are 
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therefore  convenient  to  use  In  our  problem. 

The  covariance  function  of  the  first-order  Gauss^larlcov 


process  Is  the  simple  exponential  function 

...  o  -Pi|t| 


(4.10) 


The  process  is  defined  by  Che  differential  equation 
X  +  •  w 


(4.11) 


where  w  Is  white  noise.  Because  of  its  simplicity,  this  process  is 
often  used  to  model  physical  phenomena.  It  is  not  a  consistent  model 
with  respect  to  the  gravity  field.  But  since  consistency  is  provided 
by  the  geodetic  model  this  is  not  a  major  concern.  The  essential 
parameters  of  this  process  are 


1  a  1  1  a  o  0.693 

in  2  =  ^  2  = 


(4.12a) 


(In  2)2  =  0.480 


It  should  be  noted  that  in  the  navigational  literature  the  correlation 
distance  t  is  used  instead  of  c  •  It  is  defined  as 


K(C)  =  t  Kq 


(4.13) 


where  e  is  the  basis  of  the  natural  logarithm.  In  order  to  differen¬ 
tiate  between  t  snd  t  without  too  much  verbal  acrobatics  we  will 
continue  to  call  C  the  correlation  length  while  C  will  be  called  the 
correlation  distance.  For  the  above  case 


(4.12b) 


The  covariance  function  of  the  second-order  Gauss-Markov 


process  is  given  by 


K(t)  -  a^Ml  +  e2|T|}  e 


(4.14) 
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The  corresponding  differential  equation  is  of  the  form 

X  +  +  Bj^x  -  w  (4.15) 

The  essential  parameters  are 


y  (4.16) 


J 

where  is  obtained  by  solving  the  equation 

ind  +  6242)  “  ^ 

and  C2  by  solving 


1.678 
^2  •  82 

-  2.146 

''2  -  82 

X2  -  2.817 


in(l  +  ^2^2^  *  ^2^2  “  2  . 

Similarly »  the  covariance  function  of  the  third-order  Gauss- 
Markov  process  is  given  by 

K(t)  -  03^  {1  +  B3ITI  +  i  83^  It^]}  +  m2  .  (4.17) 

with  the  corresponding  differential  equation 


X  +  363  X  +  h  ^  ^ 

and  the  essential  parameters 

S  ■<’3 

.  _  2.330 

'3  — 

r  2.903 

'3  --5— 


^3 


(4.18) 


y  (4.19) 


1.810 


J 


1 
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The  correlation  distance  t  is  computed  by  solving 

in  (1  +  C333  +  -j  C3  Sj  )  ■  ^3^3  ^  ^ 

and  the  correlation  length  C  by  solving 

In  (1  +  +  I  6^  )  -  C363  -  An  2  . 

Both  the  second  and  the  third-order  Gauss-Markov  models  are  self- 
consistent  models  for  the  representation  of  the  anomalous  gravity 
field.  Jordan  (1972)  advocates  the  third-order  model  for  the 
anomalous  potential  because  it  results  in  a  gravity  anomaly  covari¬ 
ance  function  with  a  realistic  behaviour  at  ^  »  0  or  t  »  0  respec¬ 
tively.  As  mentioned  above,  self-consistency  of  the  navigational 
model  is  not  a  major  concern  in  our  case.  The  selection  of  an 
appropriate  model  should  therefore  be  based  on  the  agreement  between 
the  essential  parameters  of  the  geodetic  model  and  the  navigational 
model.  Equations  (4.12a),  (4.16)  and  (4.19)  show  that  variance  and 
correlation  length  are  variable  while  the  curvature  parameter  is 
constant.  This  means  that  any  parameters  ,  reqxiired  by  the 

geodetic  model,  can  be  accommodated  by  the  navigational  model.  The 
selection  of  an  optimal  fit  can  therefore  be  done  on  the  basis  of 
X  only.  In  our  case  the  curvature  parameters  of  the  deflection 
covariance  functions  C(CfC)  and  C(n>n)  are  required.  They  can  be 
derived  from  equation  (4.1)  by  using  the  well  known  relations 


C 


vy  Bip 


1  3T 
VY  cos  (p  3X 


to  obtain 

cov(5p,Cp) 


1  1  3^K(P.Q) 

'p  'q  Yp  Yq  ^♦p  % 


(4.20a) 


The  differentiations  of  K(P,Q}  can  be  done  by  using 
cost^pq  *  sln4>p  s±n<^^  +  cos^^  cos^fr^  cosCX^-Xp) 

and  fonsulas  (4.4)  and  (4.3).  Then  introducing  again  planar  approxi¬ 
mation,  a  double  differentiation  of  the  functions  (4.20)  with  respect 
to  p  will  give  ic  and  thus  x  •  To  avoid  the  cumbersome  differentia¬ 
tions,  it  has  been  assumed  that  the  curvature  parameter  of  the  gravity 
anomaly  covariance  functions  is  close  enough  to  the  corresponding 
parameter  of  the  deflection  covariance  function.  Thus,  equation 
(4.8)  has  been  used.  The  actual  computation  and  comparison  of  these 
paraMters  is  given  in  section  6. 


5.  THE  STEADY  STATE  POSITION  ERROR  INDUCED  BY  GRAVITY  FIELD 
APPROXIMATIONS 

,  The  procedure  outlined  In  section  3  led  to  the  formula 

-88^^^  “  ♦2l(t,to)  ®n(t,to)  (3.16) 

for  the  steady  state  errors  of  a  linear  system.  Using  equation  (3.18) 
this  can  be  rewritten  as 

”  £2l(t,to)  iLct.to)  (5.1) 

To  obtain  t;he  steady  state  position  errors  Induced  by  different  gravity 
field  approximations,  the  system  of  differential  equations  for  this 
problem  has  to  be  set  up.  For  a  three-axes  inertial  navigation  system 
the  minimum  system  is  of  order  12  containing  the  following  elements  in 
the  state  vector 

Gj).  SPjj.  5Pg.  5Pjj.  SVjj,  6Vg,  5 ,  n,  5N)  (5.2) 

where  the  subscripts  N,  E,  D  denote  the  directions  of  the  three  axes 
of  a  local -level  system  (North,  East,  Down)  and 
c  ...  are  attitude  errors 

6p  ...  are  position  errors 

6v  ...  are  velocity  errors 

C,  n  •••  Sive  the  deflections  of  the  vertical  in  north  and  east 

direction 

6N  ...  is  the  geoidal  undulation. 

The  dynamics  matrix  for  the  first  nine  states  can  be  found  in  Britting 
(1971)  while  for  the  states  representing  the  anomalous  gravity  field 
one  of  the  models  in  section  4  can  be  used.  The  solution  of  this 
system  applying  Laplace  transform  techniques  is  a  formidable  task.  To 
give  an  idea  of  the  type  of  analytical  manipulations  necessary,  refer¬ 
ence  is  made  to  Wong  and  Schwarz  (1979)  where  the  transition  matrix 
^22  been  derived  for  the  first  9  states  in  equation  (5.2).  Con¬ 
sidering  that  the  three  additional  states  will  complicate  the  deriva¬ 
tion  considerably  and  that  the  determination  of  ^21  much  more 
complicated  than  that  of  ^22 f  decided  to  consider  only  the 


single  axis  case,  l.e«  to  completely  decouple  the  three  channels. 
Resiilts  obtained  from  numerical  Integration  techniques  Indicate  that 
this  simplified  case  represents  the  error  behaviour  of  the  more 
general  model  quite  well.  Thus  results  obtained  In  this  way  will  In 
general  be  also  valid  for  the  more  extensive  model.  The  obvious 
advantage  Is  that  a  much  smaller  system  has  to  be  dealt  with.  This 
Is  probably  the  main  reason  why  all  previous  Investigations  followed 
this  route. 

The  single  axis  case  for  one  of  the  horizontal  channels  Is 
given  by  the  differential  equations 

-K6v  -  a)^6p  -  yC 


6v 

6P 

i 

where 


G 

Ui 

where 

R 

K 

d 


6v 

-6,5  +  u, 


(5.3a) 

(5.3b) 

(5.3c) 


Is  normal  gravity  which  for  this  error  analysis  can  always 
be  replaced  by 

the  global  mean  value  of  gravity, 
is  the  Schuler  frequency 


G  V 

R 


(5.4) 


Is  the  mean  radius  of  the  earth 

is  the  feedback  gain  which  Is  related  to 

the  loop  damping  parameter  by 

K  -  2d(D  , 


(5.5) 


5  ...  Is  the  deflection  of  the  vertical  (either  5  or  n) 

and  U|  ...  Is  white  noise  defined  by 

E{uj(t)  U|(t)}  “  q^^Ct-x)  (5.6) 

where  the  power  density  spectrum  in  case  of  equation  (5.3c)  Is  given 
by 

-  20^  (5.7a) 

and  8^  and  are  defined  as  in  equation  (4.10).  The  power  density 
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spectrum  for  the  second  and  third  order  Gauss-Markov  process  are 
given  by 

92  “  (5.7b) 

and  qj  -  I6S3®  (5.7c) 

respectively.  It  should  be  noted  that  errorless  velocity  measurements 
have  been  assumed  to  avoid  a  mixing  of  error  sources  for  the  position* 
Ing  error. 

Putting  equations  (5.3)  Into  the  standard  form  (2.5)  we  have 


where  the  index  of  @  has  been  omitted  since  no  confusion  is  possible. 
Inversion  of  this  matrix  results  In 


where 

Ej  -  s^  +  sK  +  0)2  ,  (5.8b) 

Using  equations  (3.15d)  and  (3.18)  we  get 

*22(t,to)  -  l"^{(sI-F)"^}  -  (5.9) 

which  gives  us  one  of  the  required  components  In  equation  (5.1).  To 
obtain  £21  have  first  to  form  the  matrix  product 


\ 
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A  -  (sI-F)"^  £95^ 
which  Is  of  the  following  form 


gG^s^  BG^s  _  BGs 

(*2-e2)E^E^  (82-02)EjE2  (s2-62)Ej 


(82-62)EjE2  (82-62)EiE2  (s2-62)Ej 


(5.10a) 


(s2-e2)  E2 


(s2-62)  E2 


(s2-62) 


where 


j2  -  sK  +  0)2 


(5.10b) 


and  E^is  again  defined  by  equation  (5.8b).  £21  then  obtained  by 

the  inverse  Laplace  tr^sform 


£21  ^  L  ^{A} 


(5.11) 


where  A  Is  the  matrix  on  the  right-hand  side  of  equation  (5.10a).  The 
steady-state  position  error  is  determined  by  the  second  diagonal  ele¬ 
ment  in  the  P  -matrix  (5.1).  Using  the  individual  elements  of  the 
— ss 

matrices  A  and  B  defined  by  equations  (5.11)  and-  (5.9)  and  the  lin¬ 
earity  property  of  the  Laplace  transform,  we  can  write  the  variance 
of  the  position  error  In  terms  of  inverse  Laplace  transforms 


a2  -  2<Tj^  [L"Ma2i}L"^{b2i}  +  {a22>L“^ {b22} 

+  L  ^{a23}L  ^{b23)]  • 


(5.12) 


The  determination  of  the  individual  terms  in  equation  (5.12)  is  most 
simply  accomplished  by  using  the  complex  inversion  formula 


■  -  f 

2it1  J 


f(8)  ds 


(5.13 


where 


t  >  0 

8  *  u  +  iv 

i-  /T 

and  the  integration  is  performed  along  a  line  s  -  y  in  the  complex 
plane.  The  real  number  y  is  chosen  so  that  s  =  y  lies  to  the  right  of 
all  the  singularities  but  is  otherwise  arbitrary.  In  practice,  the 
integral  is  evaluated  by  considering  the  integral  enclosed  by  the 
Bromwich  contour  C 

2^  9  ds  (5.14) 


This  Integral  is  related  to  equation  (5.13)  by 


x(t)  ■  11a 


e®*^  f(s)  ds 


(5.15) 


where  R  is  the  radius  of  the  circle  and  F  the  circular  part  of  the 
contour  C.  A  solution  to  (5.14)  can  be  obtained  using  the  residue 
theorem  which  states:  If  f(s)  is  analytic  within  and  on  the  boundary 
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of  C  of  a  region  except  at  a  finite  number  of  poles  p^  within  the 


region  having  residues  p^^  respectively »  then 
f(3)  ds  ■  2Tri  Z  p^^  , 


-j  ^ 


(5.16) 


i.e.  the  solution  can  be  obtained  by  determining  the  residues  of  £(s) 
at  the  poles  py  They  can  be  found  by  the  general  formula 


S^Pj 


(5.17) 


where  n  is  the  order  of  the  pole.  For  simple  poles  eqiiation  (5.17) 
reduces  to 


p_.  »  lim  {(s-pp  f(s)}  . 

s^p .  ^ 

J 

st 

Replacing  f(s)  by  e  f(s)  and  assuming 
lim  ^2^  ^  ® 

the  solution  of  equation  (5.15)  can  be  written  as 
x(t)  -  J  P  ^ 

j 


(5.18) 


(5.19) 


(5.20) 


where  p_j  are  the  residues  of  e®'‘f(s)  at  poles  of  f(s).  Sufficient 
conditions  under  which  equation  (5.19)  is  correct  are  e.g.  given  in 
Doetsch  (1971) . 

Returning  to  equation  (5.12)  we  start  with  the  first  term 


L  { aj 1 )  *  L 


-1  r  ec^s  , 

\  (82-82)EiE2  , 

-1  f _ 

\  (s2-82)(32+g| 


ec^s _ 

^+8K+(i»2  )  (  s^-sK+u^  ) 


E^  and  can  be  factored  into 

Ej  ■  (s-fa)  (s  +  b) 

Ej  ■  (a-a)  (s-b) 


(5.21) 
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where  ______ 

a  -  -  u2 

b  -  jK  +  / -  u^' 

Because  d  <  1  in  equation  (3«5)  we  get  the  complex  roots 
a  -  |k  -  ID 

b  -  jK  +  ID 
where 


(5.22a) 


(5.22b) 


(5.22c) 


is  now  always  positive.  Using  (5.21),  we  can  write  the  inverse 
Laplace  transform  of  aai  as 


L  ^{a2i}  *  L"^  < 


_ BG^s _ 

(s^-B^) (s^-a^) (s^-b^) ^ 


Thus>  in  this  case  the  function  f(s)  has  the  six  simple  poles 
±S,  ±a,  ±b.  Using  equations  (5.20)  and  (5.18)  we  obtain 


X21<t) 


L"^{a2i} 


*  lim 
s-^  -B 

+  lim 
s-^  B 


+  lim 
s^-a 

+  lim 
s^a 


+  lim 
8^-b 


+  lim 
g-^b 


V 

(■ 


BG^s 

St 

e 

(s- 

B) (s2-a2 

)(s2 

-b2) 

BG^s 

(s+6) (s^-a^ 

’)(s2 

-b2) 

BG^s 

eSt 

(82 

-B2) (s-a 

l)  (s2 

-b2) 

BG^s 

(S2 

-B2)(s+a)(s2 

-b2) 

BG^s 

eSt 

(s2 

-b2)(s2- 

■a2)( 

s-b) 

BG^s 

(s^-B^) (s^-a^) (s+b) 


x 
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Evaluating  these  expressions  leads  to 


L"^{a2ii  “  -BG2 


[_£ 

1.(82- 


cog  h  Bt  ^  cos  h  at 
a^) (a^-6^) (a^-b^) 

cos  h  bt 


(5.23a) 


(b2-e2)(b2-a2) 

A  check  of  this  result  can  be  obtained  by  splitting  a^^  into 
ajjCs)  -  fj(s)  •  fjCs) 
where 

f^(s)  -  - - - 

(s2-a2) (s2-b2) 


6 


s2-62  ’ 


finding  the  inverse  Laplace  transforms  x^Ct)  and  x^Ct)  from  a  table 
and  applying  the  convolution  theorem 

XjiCt)  =■  x^(t)  *  XjCt) 

-t 


/ 


X^(t)  X,^(t-T)  dT 


This  is,  however,  a  more  laborious  procedure. 

The  next  term,  the  inverse  Laplace  transform  of  3^2  can  be 
obtained  by  observing  that 

L“^{a22}  *  -^321}  . 


Since  in  general 

l"^{  i  f(8)} 


r« 


(t)  dT 


where  x(t)  is  the  inverse  transform  of  f(s),  we  have  in  our  case 


L-^{a22} 


-/  ^21 


(t)  dT 
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where  is  given  by  the  right-hand  side  of  equation  (5.23a),  The 

Integrals  are  simple  to  evaluate  resulting  In 


L’^{a22> 


BG2 


sin  h  6t 


sin  h  at 


(^B(02-a2)(02_b2)  a(a2-62)(a2-b2) 


+ 


sin  h  bt  I 

b(b2-62)(b2-a2)  j 


(5.23b) 


A  check  Is  again  possible  by  use  of  the  complex  inversion  formula  and 
the  residue  theorem.  Applying  similar  techniques  for  the  other  four 
terms,  we  obtain 


L  {323} 


6G 

(a2-e2)(b2.g2) 


<(a+6)  cosh  6t 


(g2+ab)  slnh  6t  ^  (a^-B^)  e 
B  (e-b) 

(bMii^l 


L  ^  ■p  sin  Dt  , 


(5.23c) 


(5.23d) 


compare  equation  (2.15c), 


L“^{b22> 


-JjKt 


and 


L-'{b„) 


-Ci 


K  sin  Dt 
2D 


-at 


+  cos  Dt 


-bt 


(b-a)(6-^)  (a-b)(B-b)  (a-B)(b 


-Bt  1 

a 

= -  k 

i)(b-B)^ 


(5.23e) 


(5.23f) 


The  variance  of  the  position  error  can  now  be  obtained  by  performing 
the  multiplications  in  formula  (5.12)  using  the  expressions  (5.23a)  to 
(5.23f).  The  somewhat  lengthy  formula  is  not  written  out  here.  It 
gives  the  steady  state  error  as  a  function  of  time.  Often  we  are  only 
interested  in  the  behaviour  of  the  error  after  it  has  settled  down, 
i.e.  for  t  -►  »  .  The  relevant  expressions  have  been  published  by 


several  authors,  see  e.g.  Jordan  (1973),  and  are  not  rederlved  here. 
They  are  usually  obtained  by  applying  the  final  value  theorem  of 
Laplace  transforms  directly  to  equation  (5.12). 

For  the  first-order  Gauss-Markov  model  with  the  covariance 
function  (4.10)  we  obtain 


K  +  B, 


(B^^+KBi+  0)2) 


(5.24) 


For  the  second-order  Gauss-Markov  model  with  the  covariance  function 


(4.14)  we  get 


P2  2 


282^+62)^  +  Ka)2 
(82^  +  K82  + 


(5.25) 


Finally,  use  of  the  third-order  Gauss-Markov  model  with  covariance 
function  (4.17)  results  in 


2-2 
P3  3 


'(7KB^^  +3^k2  +Ka)2  +5B3^)  +  63^(83+ K)^ 
(83^  +  KB3  +  to2)^ 


(5.26) 


Equations  (5.24)  to  (5.26)  show  that  the  elimination  of  time  from 
the  error  function  results  in  a  considerable  simplification  of  the 
expressions.  When  comparing  the  effect  of  different  gravity  field 
approximations  we  are  primarily  interested  in  the  average  behaviour 
of  the  error  over  long  time  intervals.  Thus,  equations  (5.24)  to 
(5.26)  will  serve  our  purpose  well  and  will  therefore  be  used  in  the 
next  section* 

The  discussion  so  far  has  been  restricted  to  the  two  hori¬ 
zontal  channels,  i.e.  only  the  position  errors  due  to  unmodelled 
deflections  of  the  vertical  have  been  considered.  Formulas  for  the 
height  error  due  to  unmodelled  gravity  anomalies  or  geoidal  undula¬ 
tions  can  be  derived  in  an  analogous  way.  They  can  for  instance  be 
found  in  Jordan  (1973).  Due  to  the  fact  that  the  natural  frequency  of 
the  vertical  channel  is  higher,  these  errors  are  much  smaller,  usually 
less  than  5%  of  those  generated  in  the  horizontal  channels.  They  are 
therefore  not  given  here. 


The  same  procedure  which  was  used  to  derive  the  steady  state 
position  error  can  be  employed  to  determine  the  steady  state  velocity 
error  induced  by  insufficient  deflection  information.  In  this  case, 
the  first  diagonal  term  of  the  P^^'-matrix  contains  all  the  necessary 
information.  Formulas  can  be  found  in  Levine  and  Gelb  (1969)  and 
Jordan  (1973) . 
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6.  RESULTS 

To  use  the  formulas  developed  in  the  last  section,  three 
steps  are  necessary.  First,  the  different  gravity  field  approxi¬ 
mations  have  to  be  specified  using  the  geodetic  model  described  in 
section  4.  Second,  the  essential  parameters  of  the  resulting  co- 
variance  functions  have  to  be  determined.  Finally,  a  suitable  navi¬ 
gational  covariance  model  must  be  selected  by  fitting  the  essential 
parameters  of  the  geodetic  model  to  those  of  the  navigational  model. 
After  that,  formulas  (5.24)  to  (5.26)  can  be  used  to  determine  the 
position  error  induced  by  insufficient  knowledge  of  the  deflections 
of  the  vertical. 

Gravity  field  approximations  of  increasing  accuracy  can  be 
represented  by  spherical  harmonic  expansions  of  increasing  degree 
and  order  N.  When  using  such  models  in  an  airborne  inertial  navigator, 
the  long  wave-length  features  of  the  gravity  field  are  eliminated. 

This  means  that  with  Increasing  M  the  unmodelled  effects  become  more 
and  more  local  in  character.  The  covariance  representation  of  geo¬ 
potential  models  of  arbitrary  degree  and  order  N  can  be  simulated  by 
subtracting  the  first  N  coefficients  from  the  model  (4.1).  Covari¬ 
ance  computations  with  such  models  are  done  in  exactly  the  same  way 
as  with  the  full  model.  We  thus  have  a  consistent  and  simple  way  to 
characterize  different  gravity  field  approximations. 

Table  6.1  lists  the  models  used  in  the  following  computations 
and  indicates  reasons  for  their  choice.  The  first  three  models 
represent  gravity  field  approximations  available  today.  The  next 
three  will  hopefully  become  available  on  a  global  scale  in  the  not  too 
distant  future.  The  last  two  represent  models  of  well  surveyed  areas 
but  cannot  be  considered  as  given  globally.  Since  gravity  field  infor 
mation  is  often  given  in  gridded  form,  the  grid  interval  corresponding 
to  a  geopotential  model  of  degree  and  order  N  is  also  listed  in  the 
last  column.  It  has  been  computed  by  the  approximate  formula 

grid  interval  »  ir/N  , 


which  is  useful  for  discrete  function  values  but  has  limitations  when 


applied  to  mean  values,  see  e.g.  Rapp  (1977). 


Degree  and  order  N 
of  geopotential  model 

Significance 

gridded 

equivalent 

N  -  2 

Normal  gravity  model 

Geopotential  model  derived 

N  •  18 

from  orbital  analysis  of 

10“ 

satellites. 

N  »  36 

Geopotential  model  from  a 
combination  of  orbital 
analysis  of  satellites  and 
terrestrial  gravimetry. 

B 

N  -  90 

Geopotential  model  from  a 
combination  of  satellite 
altimetry  and  terrestrial 
gravimetry. 

2“ 

N  -  180 

Geopotential  model  from  a 
combination  of  terrestrial 
gravimetry  and  satellite  to 
satellite  tracking  or  satel¬ 
lite  gradiometry. 

1® 

N  »  300 

as  N  =  180,  using  more 
optimistic  assumptions . 

36' 

N  »  720 
N  -  1080 


Local  gravity  field  repre¬ 
sentations  obtained  from 
terrestrial  gravimetry  or 
airborne  gradiometry. 


Table  6.1  :  Gravity  Field  Approximations  Used 


The  gravity  field  models  liste<^  in  table  (6.1)  are  needed  at 
flying  altitudes.  As  a  typical  range,  flying  heights  between  10  km 
and  20  km  have  been  chosen.  The  upward  continuation  of  the  covariance 
model  (4.1)  is  obtained  by  an  appropriate  change  of  the  radial  dis¬ 
tances  rp  and  •  Again,  the  covariance  computation  is  performed  by 
the  same  algorithm  as  before. 
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The  basic  covariance  model  used  to  compute  the  models  in 
table  6.1  has  the  p  .ameters 

A  ■  607.57  mgal^ 

k  -  - - - -  (6.1) 

(n-l)(n-2)(n+24) 

s  -  0.998444 

where 


and  ^  *  6371  km 

is  the  mean  radius  of  the  earth.  As  usual,  the  parameter  A  corres¬ 
ponds  to  the  gravity  anomaly  covariance  function.  Its  essential 
parameters  are 

C  •  1762  mgal^ 
o  ® 

t  -  4814 
X  -  13.7 

They  correspond  somewhat  better  to  the  observational  data  than  the 
parameters  given  by  Tscheming  and  Rapp  (1974).  For  a  discussion, 
see  Schwarz  and  Lachapelle  (1980) . 

The  different  models  In  table  6.1  were  computed  using  the 
subroutines  COVAX  by  Tscheming  (1976)  and  COVAPP  by  Siinkel  (1979). 

A  graphic  representation  of  the  results  is  given  in  figures  6.1  to 
6.4.  They  show  the  covariance  fxmctions  for  the  two  deflection  com¬ 
ponents  at  altitudes  h  »  10  km  and  h  «■  20  km.  The  gravity  field 
approximation  is  characterized  by  the  degree  and  order  N  of  the  geo¬ 
potential  model.  The  covariance  functions  N  >  180  for  h  «  10  km  and 
N  >  90  for  h  »  20  km  are  too  small  to  be  represented  in  the  given 
scale. 

Three  distinct  features  can  be  observed  from  these  graphs. 
First  of  all,  the  covariance  function  of  the  ^'-component  is  much 
smoother  than  that  of  the  C^component  (note  that  the  scales  of  the 


V. 


(arc  aacV 
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eOVABI*NCE[i,f]  »T  »-’<»“ 


COVARIAMCE[n,.]  at  »=10<»" 


Fig.  6*2 


Deflection  covariance 


function  C(n,n) 


h  •  10  to 


arc  min(^) 


Fig*  6*4  Deflection  covariance  function  C(ri»ri)  at  h  •  20  km 
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axis  are  different  between  the  two  seta  of  figures) .  This  means 
that  we  can  expect  a  major  difference  in  the  correlation  length. 

Second,  the  effect  of  increasing  the  degree  and  order  N  of  the  approxi* 
mation  is  a  reduction  of  both  the  variance  and  the  correlation  length 
of  -the  covariance  function.  This  feature  is  easily  explained  because 
the  gravity  field  represented  by  the  covariance  function  becomes  more 
and  more  local  with  increasing  N.  Finally,  the  effect  of  upward 
contintiatlon  is  a  decrease  in  the  variance  but  an  increase  in  the 
correlation  length.  Again,  this  is  physically  meaningful  because  the 
high  frequency  spectrum  is  more  affected  by  the  attenuation  with 
height.  It  should  be  noted  that  these  results  cannot  be  generalized 
to  other  functionals  of  the  anomalous  potential,  as  e.g.  the  geoidal 
undulations  or  second-order  gradients.  Since  they  are  influenced  by 
different  parts  of  the  spectrum  their  characteristical  behaviour  with 
respect  to  changes  in  N  and  h  is  quite  different.  For  a  discussion 
of  this  point,  see  Schwarz  (1981). 

The  second  step  to  the  use  of  equations  (5.24)  to  (5.26)  is 
the  determination  of  the  essential  parameters  for  all  covariance 
functions.  They  can  be  obtained  by  applying  the  formulas  of  section  4 
to  the  output  of  the  subroutines  mentioned  above.  Results  are  summar¬ 
ized  in  tables  6.2  to  6.5  which  correspond  to  figures  6,1  to  6.4. 

Besides  the  features  already  mentioned,  the  change  in  the 
curvature  parameters  x  is  of  specific  interest.  It  decreases  with 
increasing  N  and  increasing  h.  Since  the  decrease  with  N  is  especially 
pronounced  for  the  first  few  degrees,  all  models  with  N  >  18  are  quite 
well  represented  by  covariance  functions  of  second  and  third  order 
Gauss-Markov  models.  This  is  an  important  result  because  it  secures 
a  good  fit  of  the  geodetic  to  the  navigational  models.  It  should  be 
noted,  however,  that  this  is  not  true  any  more  for  inertial  navigation 
on  the  surface  of  the  earth.  Due  to  the  increase  of  x  with  decreasing 
h,  the  curvature  parameters  tend  to  be  too  large  for  low-order  Gauss- 
Markov  models.  That  this  is  not  a  feature  of  the  specific  geodetic 
covariance  model  used  can  be  gathered  from  a  table  of  empirically 

1 
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Model 

N  - 

Variance 

(arc  sec)^ 

CorreJ 
length  c 
(arc  min) 

.ation  _ 

distance  Z 
(arc  min) 

Curvature 

parameter 

X 

2 

32.5 

81.2 

136. 

7.41 

18 

19.8 

40.3 

54.2 

3.94 

36 

14.9 

30.5 

38.9 

3.00 

90 

7.94 

19.3 

23.3 

1.87 

180 

3.57 

12.7 

14.9 

1.48 

300 

1.47 

9.03 

1 

10.5 

1.32 

720 

0.106 

4.61 

5.29 

1.20 

1080 

1 

1 

0.0137 

3.27 

3.75 

1.23 

Table  6.2  :  Essential  Parameters  of  Deflection 

Covariance  Function  C(?»5)  for  h  *  10  km 


Model 

N  - 

Variance 

(arc  sec)^ 

Correl 
length  c 
(arc  min) 

ation 

distance  C 
(arc  min) 

Curvature 

parameter 

X 

2 

32.5 

211. 

387. 

7.41 

18 

19.8 

92.0 

132. 

3.94 

36 

14.9 

66.4 

89.9 

3.00 

90 

7.94 

39.2 

49.7 

1.87 

180 

3.57 

24.7 

30.2 

1.48 

300 

1.47 

17.1 

20.5 

1.32 

720 

0.106 

8.51 

10.0 

1.20 

1080 

0.0137 

i 

_ 1 

6.00 

7.05 

1.23 

Table  6.3  :  Essential  Parameters  of  Deflection 

Covariance  Function  CCn^n)  for  h  10  km 
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Model 

N  - 

Variance 

(arc  sec)^ 

Cor  re. 
length  c 
(arc  min) 

Lation 

distance  t 
(arc  min) 

Curvature 

parameter 

X 

2 

113. 

182. 

6.50 

18 

52.5 

67.8 

3.02 

36 

38.8 

47.9 

2.22 

90 

4.35 

23.5 

27.8 

1.57 

180 

1.42 

15.0 

17.3 

1.32 

300 

0.390 

10.3 

11.8 

1.24 

720 

0.00714 

4.94 

5.69 

1.17 

1080 

0.000300 

3.50 

4.00 

1.16 

Table  6.4  :  Essential  Parameters  of  Deflection 

Covariance  Function  C(^,^)  for  h  ^  20  km 


Model 

Variance 

Correlation 

— 

Curvature 

length  c 

distance  c 

parameter 

N  - 

(arc  sec)^ 

(arc  min) 

(arc  min) 

X 

2 

290. 

511. 

6.50 

18 

116. 

159. 

3.02 

36 

81.3 

106. 

2.22 

90 

46.2 

57.1 

1.57 

180 

28.3 

34.0 

1.32 

300 

0.390 

19.1 

22.6 

1.24 

720 

0.00714 

9.14 

10.7 

1.17 

1080 

0.000300 

6.30 

7.50 

1.16 

Table  6.5  :  Essential  Parameters  of  Deflection 

Covariance  Function  C(n»n)  for  h  -  20  km 
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determined  x^values  published  in  Schwarz  and  Lachapelle  (1980).  They 
have  an  average  alze  of  x  *  ^  when  referenced  to  a  geopotential  model 
of  degree  and  order  22.  This  means  that  at  h  »  0  the  simple  naviga¬ 
tional  models,  which  typically  have  X'^values  below  3,  represent  only 
local  disturbances  of  the  anomalous  gravity  field  well.  They  will 
therefore  give  reliable  estimates  only  when  used  with  a  gravity  field 
approximation  of  high  degree  and  order.  Since,  at  present,  geopoten¬ 
tial  models  with  a  large  enough  N  are  not  available,  navigational 
models  are  needed  which  admit  a  larger  curvature  parameter.  It  is 
also  interesting  to  note  that  for  the  same  reason  the  covariance 
function  of  the  first-order  Gauss-Harkov  process  is  not  a  good  model 
of  the  anomalous  gravity  field  anywhere  in  the  range  0  1  h  £  20  km. 

The  third  step  to  the  use  of  equations  (5.24)  to  (3.26)  is  the 
fitting  of  the  geodetic  models  to  one  of  the  navigational  models  by 
way  of  the  essential  parameters.  This  amounts  to  obtaining  a  good  fit 
for  the  curvature  parameter  because  the  variance  and  the  correlation 
length  of  the  two  models  can  always  be  made  to  agree.  Using  the  geo¬ 
detic  model  parameters  given  in  tables  (6.2)  to  (6.3)  a  second  or 
third  order  Gauss-Markov  process  is  chosen  depending  on  whether  the 
X-parameter  in  formula  (4.16)  or  (4.19)  is  closer.  Then,  transforming 
C  into  the  time  domain  by 


where  v  is  the  aircraft  velocity,  formula  (5.25)  or  (5.26)  can  directly 
be  applied.  For  a  definition  of  terms  in  these  equations,  see  formulas 
(5.3)  to  (5.5).  A  damping  factor  of  0.3  has  been  used  throughout. 

Results  for  aircraft  velocities  of  300  knots,  500  knots,  and 
800  knots  are  presented  in  tables  6.6  and  6.7.  They  show  clearly  that 
the  single  most  influential  factor  for  a  reduction  of  position  errors 
is  the  degree  of  the  gravity  field  approximation  model.  Flying  alti¬ 
tude,  aircraft  velocity,  and  correlation  distance  have  some  effect  but 
are  not  decisive.  Using  one  of  the  available  geopotential  models  of 
degree  and  order  36  instead  of  the  normal  gravity  field  would  reduce 
the  errors  by  about  40%  on  average,  and  by  50%  in  case  of  a  high  speed 


Steady-state  position  error  induced  by 


C(€,C) 

Model  300  kn  500  kn  800  ka 

N  • 


197 

208 

207 

162 

149 

128 

135 

117 

98 

69 

53 

41 

37 

28 

22 

19 

15 

12 

4 

3 

3 

2 

1 

1 

C(n,n) 

300  kn  500  ka  800  ka 


C(C,S)  and  C(n,n) 
300  kn  500  kn  800  kn 
m) 


267 

279' 

284 

224 

220 

206 

194 

183 

131 

119 

94 

73 

65 

50 

39 

34 

23 

21 

6 

5 

4 

3 

2 

1 

Table  6.6  :  Gravity  Induced  Position  Error  for  h  =«  10  km 


Table  6.7  :  Gravity  Induced  Position  Error  for  h  «  20  km 


Model 
N  - 

300  ka 
(m) 

Steady-stai 

C(5,C) 

500  kn  800  kn 
(m)  (m) 

te  position  error  a  i 
P 

c(n,n) 

300  kn  500  kn  800  kn 
(m)  (m)  (m) 

nduced  by 

C(c,5)  and  C(Ti,n) 

300  kn  500  kn  800  kn 
(m)  (m)  (m) 

2 

174 

184 

191 

163 

166 

173 

238 

248 

258 

18 

140 

134 

119 

130 

137 

140 

191 

192 

284 

36 

108 

87 

68 

119 

116 

101 

161 

145 

122 

90 

56 

43 

33 

74 

62 

49 

93 

75 

59 

180 

25 

19 

15 

35 

27 

21 

43 

33 

26 

300 

14 

8 

6 

15 

11 

9 

21 

14 

11 

720 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1080 

0.2 

0.2 

0.1 

0.3 

0.2 

0.2 

0.4 

0. 

3  0.1 

1 

aircraft.  With  gravity  models  expected  to  become  available  in  the 
near  future  (N  “  180) ,  the  position  error  could  be  reduced  to  about 
30  m  to  50  m,  i.e.  to  13%  or  20%  of  what  it  is  today.  From  there, 
progress  will  be  slow.  To  reach  the  meter  range,  a  gravity  field 
approximation  equivalent  to  an  expansion  of  degree  and  order  1000 
would  be  necessary.  This  result  is  not  surprising.  It  demonstrates 
the  well-known  fact  that  the  medium  and  high  frequency  spectrum  con¬ 
tributes  considerably  to  the  deflections  of  the  vertical.  Or,  in 
other  wo;rds,  the  relative  contribution  of  local  effects  is  not 
negligible. 

Aircraft  speed  affects  the  position  error  to  some  extent. 

For  all  but  the  normal  gravity  field  approximation  a  high  aircraft 
velocity  will  improve  results.  The  reason  for  this  can  be  found  in 
the  dependence  of  the  position  error  on  6.  The  error  curve  increases 
to  a  6  -  value  of  about  3h  ^  or  4h  ^  and  then  falls  off  slowly,  see  e.g. 
figure  4  in  Jordan  (1973).  Due  to  the  long  correlation  distances  of 
the  normal  model,  a  3  -  value  of  3.5h  and  thus  the  largest  position 
error,  is  reached  at  about  v  640  knots  for  5  and  v  =  1790  knots  for 
n*  For  all  the  other  models  a  speed  of  300  knots  corresponds  to  a  3 - 
value  well  above  3.5h  Thus,  an  increase  in  v  will  result  in  a 
decrease  of  the  position  error.  It  also  explains  why  a  large  correla¬ 
tion  distance  reduces  the  error  in  the  normal  model  and  increases  it 
in  all  the  other  models.  An  increase  in  flying  altitude  will  in 
general  reduce  the  position  error.  In  this  case,  the  smaller  variance 
of  the  covariance  curve  more  than  compensates  for  the  increase  of  the 
position  error  due  to  the  larger  correlation  distance. 

The  position  errors  listed  in  tables  6.6  and  6.7  correspond 
to  an  average  global  behaviour  of  the  anomalous  gravity  field.  In 
areas  where  considerable  deviations  from  this  average  behaviour  occur, 
larger  or  smaller  errors  can  be  expected.  An  indication  of  the  possible 
range  of  values  is  given  in  table  1  of  Bernstein  and  Hess  (1976). 
Although  results  are  not  directly  comparable  because  mean  gravity  values 
and  a  somewhat  different  technique  were  used  there,  their  'best'  and 
'worst'  case  is  quite  close  to  what  can  be  expected  when  using  the 
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empirical  covariance  functions  In  Schwarz  and  Lachapelle  (1980). 

The  results  also  permit  some  conclusions  on  the  effect  of 
measuring  errors  in  the  data.  Random  errors  will  be  strongly 
attenuated  with  height  and  will  have  a  negligible  effect  on  position 
accuracy.  However,  data  errors  which  generate  spatial  correlations 
In  the  gravity  field  model  will  result  In  position  errors  of  consider¬ 
able  size.  Errors  of  this  type  are  e.g,  data  biases  over  an  extended 
area  or  blunders  In  discrete  points  which  by  way  of  upward  continua¬ 
tion  cause  correlations.  To  give  an  indication  of  the  possible 
effects,  an  error  covariance  function  with  a  variance  of  (1  arc  sec)^ 
and  different  correlation  distances  has  been  used  to  compute  the 
position  errors  In  table  6.8.  Although  they  are  at  present  smaller 


Table  6.8  :  Position  Error  Due  to  Data  Bias 


than  the  effects  coming  from  an  Insufficient  gravity  field  approxima¬ 
tion,  they  are  large  enough  to  warrant  a  careful  analysis  of  the 
existing  gravity  field  data  to  detect  errors  of  this  kind. 
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7.  C<»rCLUSI0NS  AND  REC(»1MENDATIONS 

Simplified  analytical  error  models  can  be  used  with  advantage 
to  study  the  effect  of  gravity  field  approximations  on  the  position 
accuracy  of  airborne  inertial  navigation.  In  this  approach  the  com¬ 
ponents  of  the  gravity  vector  have  to  be  expressed  as  elements  of  a 
state  vector  which  means  that  they  are  usually  modelled  as  time 
correlated  noise  processes.  Most  of  the  processes  proposed  for  this 
purpose  in  the  navigational  literature  are  of  the  Gauss-Markov  type 
and  imply  a  certain  covariance  structure  of  the  anomalous  gravity 
field.  Since  global  geodetic  covariance  models  cannot  be  cast  into 
a  form  suitable  for  the  state  space  approach,  the  question  arises 
how  well  the  navigational  models  approximate  the  geodetic  covariance 
models.  A  direct  comparison  of  the  functions  is  not  possible. 
However,  using  the  essential  parameter  approach  proposed  by  Mortiz, 
a  meaningful  comparison  can  be  made.  Results  presented  in  this 
report  show  that  none  of  the  proposed  navigational  models  gives  a 
satisfactory  fit  when  the  normal  gravity  field  is  used  as  the  basic 
approximation.  Replacing  the  normal  field  by  one  of  the  available 
geopotential  models  (IS  ^  N  36)  results  in  a  good  agreement  between 
the  geodetic  covariance  function  at  flying  altitudes  and  the  covari¬ 
ance  function  implied  by  a  third-order  Gauss-Markov  process.  For 
higher  order  reference  fields  second-order  Markov  processes  give  a 
better  fit.  The  first-order  Gauss-Markov  process,  often  used  in 
applications  because  of  its  simplicity,  is  not  a  good  model  for  the 
anomalous  gravity  field.  At  the  surface  of  the  earth,  none  of  the 
navigational  models  give  a  satisfactory  agreement  with  all  three 
essential  parameters  determined  from  actual  data.  Thus,  to  use 
these  techniques  in  the  marine  and  land-based  environment,  state 
space  models  are  needed  which  give  a  better  fit  to  empirical  covari¬ 
ance  functions. 

Use  of  the  analytical  models  show  that  the  position  error  is 
mainly  due  to  poorly  modelled  deflection  components.  The  gravity 
anomaly  and  the  geoidal  undulations  have  a  much  smaller  effect  and 


can  in  general  be  neglected.  The  size  of  the  error  is  strongly 
dependent  on  the  degree  and  order  of  the  gravity  field  model  used* 
Flying  altitude,  aircraft  velocity,  and  correlation  distance  have 
some  effect  but  are  not  decisive.  Using  one  of  the  available  geo*- 
potential  models  of  degree  and  order  36  instead  of  the  normal  gravity 
field  will  reduce  the  present  position  error  of  about  a  •  300  m  by  40% 
on  average  and  by  50%  in  case  of  a  high  speed  aircraft.  With  gravity 
models  expected  to  become  available  in  the  near  future  (N  *  180) ,  the 
position  error  could  be  reduced  to  about  15%  of  its  present  size. 

From  there  on,  progress  will  be  slow.  An  approximation  equivalent  to 
an  expansion  of  degree  and  order  1000  will  be  needed  to  reach  the 
meter  range. 

The  Implementation  of  high-order  reference  fields  in  real  time 
poses  considerable  problems  of  which  the  adequate  representation  of 
the  available  Information  is  the  most  pressing  one.  Obviously,  only 
interpolation  methods  with  local  support  at  flying  altitude  are  fea¬ 
sible  for  this  purpose.  They  will  require  storage  capacity  but  will 
not  be  demanding  in  terms  of  computing  power.  Investigations  on 
suitable  representations  of  the  gravity  vector  at  flying  altitude  and 
on  simple  Interpolation  methods  which  can  be  executed  in  real  time 
are  needed  to  make  practical  use  of  these  results. 

Finally,  errors  in  the  gravity  field  data  which  generate 
spatial  correlations  in  the  model  will  result  in  position  errors  of 
considerable  size.  Typical  errors  of  this  type  are  biases  between 
different  data  bases  and  every  effort  should  be  made  to  detect  and 
eliminate  them*  Similarly,  blunders  in  discrete  points  will  generate 
correlations  at  flying  altitude.  In  this  case,  the  method  used  for 
upward  continuation  has  some  effect  on  the  extent  of  the  correlation 
and  methods  which  reduce  the  influence  of  such  errors  should  be 
developed. 
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